setwd("/Users/xichen/Desktop/data")
library(geepack)
ds <- read.table("leprosy.txt",header=T)
#reshape data from short to long format
dslong <- reshape(ds, idvar="id", varying=c("y1","y2"), 
         v.names="y", timevar="time", time=0:1, direction="long")
dslong <- dslong[order(dslong$id, dslong$time),]
attach(dslong) 
summary(dslong)
#create interaction terms
timeA <- time*(drug=='A')
timeB <- time*(drug=='B')
#first model treat treatment A and B seperately
model1=geeglm(y ~ time + timeA + timeB, data=dslong, id = id, 
               waves=time, family=poisson("log"), corstr="exch", 
               std.err="san.se")
summary(model1)
#calculate the 95% CI and ratio
cc <- coef(summary(model1))
citab <- with(as.data.frame(cc),
              cbind(mean=Estimate,
                    lwr=Estimate-1.96*Std.err,
                    upr=Estimate+1.96*Std.err))
rownames(citab) <- rownames(cc)
ci1=exp(citab)
ci1

#second model combine treatment A and B together
timeAB <- time*I(drug != 'C')
model2=geeglm(y ~ time + timeAB, data=dslong, id = id, 
               waves=time, family=poisson("log"), corstr="exch", std.err="san.se")
summary(model2)
#calculate the 95% CI and ratio
cc <- coef(summary(model2))
citab <- with(as.data.frame(cc),
              cbind(mean=Estimate,
                    lwr=Estimate-1.96*Std.err,
                    upr=Estimate+1.96*Std.err))
rownames(citab) <- rownames(cc)
ci2=exp(citab)
ci2

detach(dslong)